SuperdifFusion in a Honeycomb Billiard 
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We investigate particle transport in the honeycomb billiard that consists of connected channels 
placed on the edges of a honeycomb structure. The spreading of particles is superdiffusive due to 
the existence of ballistic trajectories which we term perfect paths. Simulations give a time expo- 
nent of 1.72 for the mean square displacement and a starlike, i.e., anisotropic particle distribution. 
We present an analytical treatment based on the formalism of continuous-time random walks and 
explain both the time exponent and the anisotropic distribution. In billiards with randomly dis- 
tributed channels, conventional diffusion is always observed in the long-time limit, although for small 
disorder transient superdiffusional behavior exists. Our simulation results are again supported by 
an analytical analysis. 
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I. INTRODUCTION 

Billiard systems, i.e., point particles moving freely in 
areas bounded by closed curves from which they reflect 
specularly, are a paradigm of classical mechanics illus- 
trating the difference between regular and chaotic motion 
PIQ- Moreover, they are helpful to explore the relation 
between classical and quantum mechanics Q, where in- 
terference of waves becomes important as also observed 
in optical resonators |j] or by chaotic scattering in optical 
billiards [jj. 

In this paper we investigate a special example of the 
so-called infinite domain or extended billiard, where par- 
ticles' motion is unrestricted. The most famous example 
invented by Lorentz [(| to model electrons in a metal is 
the Lorentz gas, where particles reflect specularly from 
randomly distributed spherical scatterers. A periodic 
version, the so-called Sinai billiard 0, is illustrated in 
Fig. nji). For special paths, the particles possess an 
infinite horizon, i.e., they move straight without being 
scattered [see Fig. |TJi)]. These paths are responsible for 
the observation that a collection of particles with arbi- 
trary initial direction experience superdifFusion, i.e., in 
the Sinai billiard their mean-square displacement grows 
as tint , where t denotes time [8| . If sufficiently large 
scatterers are placed on a hexagonal lattice see Fig. 
|TJ>)] or if small scatterers fill the interstitial space of the 
Sinai billiard in Fig. f^i) > particles always have a fi- 
nite horizon and their motion is purely diffusive. 

Recently, transport properties in one-dimensional ex- 
tended billiards have been studied fiol JllJ . also with spe- 
cial emphasis on heat conduction |l2l Il3| | . One exam- 
ple, a periodic arrangement of channels (see Ref. [Tl|p. 
is pictured in Fig. Although all particles in such a 
channel billiard have a finite horizon, there exist paths 
illustrated in Fig. f^;) and termed perfect paths in the 
following, where the particles always move ballistically 
in one direction and therefore cause superdifFusion. 

In this paper, we study a two-dimensional extended 
billiard, where channels are placed on the edges of the 
honeycomb structure, and denote it honeycomb billiard 




FIG. 1: a) Ballistic paths in a Sinai billiard, where parti- 
cles experience an infinite horizon, give rise to superdiffusion. 
b) If the scatterers on hexagonal lattice points are sufficiently 
large, the particle trajectories always have a finite horizon and 
the spreading of particles is diffusive, c) A one-dimensional 
channel billiard with a perfect path giving rise to superdiffu- 
sion. 



(see Fig. [3| ■ The spreading of particles in such a billiard 
is also superdiffusive due to the existence of numerous 
ballistic or perfect paths examples of which are illustrated 
in Fig. |21 Note that path (3) is equivalent to the one 
in the one-dimensional billiard of Fig. f^). Our system 
is an example where particles perform a Levy walk [Trl 
llH ITsf . Very long effective steps along almost perfect 
paths lead to superdiffusion. We study it with the help 
of computer simulations and motivate the time exponent 
for the particles' mean square displacement within the 
velocity model [l9| of continuous-time random walks [TtI 
ITsj . We furthermore look at random distortions of the 
honeycomb billiard and show that for small distortions 
a transient superdiffusive regime exists whereas for large 
times and also for large distortions the spreading of the 
particles is always diffusive. Our numerical results are 
again supported by an analytical analysis. 

Originally, our work was motivated by the observation 
of photon channelling in foams 0, . Light in foams 
is reflected at the liquid-gas interface of the thin films 
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FIG. 2: Honeycomb billiard with three examples of perfect 
paths. 



which ultimately leads to a diffusive transport of photons 
through the system. Experiments and theoretical consid- 
erations show that the photons have a higher probability 
to move in the liquid phase of the films, a phenomena 
that was then termed photon channelling. The channel 
billiards studied here are an extreme case where the pho- 
tons always move in the liquid phase. 

Finally, we add a note concerning the classification of 
our billiard system. Arnol'd's famous theorem states that 
the phase space of integrable systems in classical me- 
chanics is a torus |2j|. On the other hand, the Lorentz 
gas is chaotic and therefore nonintegrable Besides 
quasi-integrable systems 2], Richens and Berry identify 
pseudointegrable systems with chaotic properties whose 
phase space is a multi-handled sphere instead of a torus 
|21| . As an example, they investigate a system similar to 
the honeycomb billiard but with the hexagons replaced 
by squares and show that the corresponding phase space 
is a five-handled sphere [2l|. Performing an analogous 
investigation for the honeycomb billiard, we find a ten- 
handled sphere as phase space |22|. 

In the following, we first introduce details of our bil- 
liard system and the method of simulations. In Section 
IIIII we report our numerical results, discuss analytical 
approaches in Section Hvl and then end with conclusions. 
Appendix contains details of the Levy- walk model for 
the honeycomb billiard. 



II. MODEL SYSTEM AND METHOD OF 
SIMULATIONS 

The objective of this paper is to study the dynam- 
ics of particles in two-dimensional channel billiards. The 
construction of the regular honeycomb billiard, as illus- 
trated in Fig. [3 is obvious. However, we also want to 
investigate random channel billiards. To create them, we 
employ Voronoi tessellations of the plane p3. |25| . They 
are generated from a distribution of seed points in a sim- 
ulation box, for which Voronoi polygons are constructed 
in complete analogy to the Wigner-Seitz cells of period- 



ically arranged lattice sites. For example, a triangular 
lattice of seed points gives the regular honeycomb struc- 
ture whose edges we choose to have the length lo. In the 
following all lengths are given in units of Iq. Then we sys- 
tematically introduce disorder by shifting the seed points 
along a randomly chosen displacement vector whose mag- 
nitude is equally distributed in the intervall [0, Sr]. All 
of our Voronoi tessellations are produced by the soft- 
ware Triangle [26j: examples are presented, e.g., in Ref. 
[27|. Typically, they contain approximately 15000 cells 
which corresponds to a quadratic simulation box with 
edge length 200 lo- This simulation box is extended in 
all spatial directions by periodic boundary conditions. 

Now, we arrive at a random channel billiard by placing 
a channel of width d on each edge of the Voronoi tessel- 
lation. Only modest disorder quantified by Sr < 0.3 is 
investigated so that all cells still have six edges. This 
avoids the situation that four instead of three channels 
meet when we construct the billiard system, which sim- 
plifies the determination of the particle path. Particles 
perform a ballistic motion with a constant velocity c in- 
side the channels; when they hit the boundary they are 
reflected specularly. In the following, we use the time 
scale lg/c to rescale time. 

Typically, we launch 10000 particles at one vertex of 
the underlying Voronoi tcsselation in an angular range of 
60° and let them run during a time t — 10 5 . At several 
times, we calculate the mean square displacement (r 2 ), 
where r denotes the position vector of the particles in 
the particle cloud, and plot it as a function of t. When 
applicable, diffusion constants in units of Iqc are then 
determined from a fit to (r 2 ) = ADt. 



III. SUPERDIFFUSION: RESULTS FROM 
SIMULATIONS 

In this section we present our numerical results for 
the exact honeycomb billiard and the random channel 
billiard. Figure |3 where we plot the mean-square dis- 
placement as a function of time for different Sr, summa- 
rizes our main results. In the exact honeycomb billiard 
(Sr = 0, plus symbols), the particles exhibit superdiffu- 
sion, i.e., in (r 2 ) oc t v the time exponent is larger than one 
and assumes the value v = 1.72 ± 0.02. Within the nu- 
merical error, the exponent is independent of the channel 
width d as long as d is small enough so that the particle's 
horizon is finite. Nevertheless, even for a finite horizon, 
we find ballistic trajectories, termed perfect paths in the 
following, in the sense that the particles move, on aver- 
age, in one direction although they experience numerous 
reflections in the channels. Examples of such trajecto- 
ries are illustrated in Fig. [2 path (3) is equivalent to the 
"propagating periodic orbit" in the work of Sanders and 
Larralde [see Fig. and Ref. P3l- 

In the framework 

of Levy walks, they can be considered as steps of infinite 
length and therefore are responsible for the superdiffu- 
sive behavior. We will investigate them in more detail 
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perfect honeycomb billard 
small disorder 8r=10" 6 
random channels 5r=0.3 




FIG. 3: Mean-square displacement as a function of time for 
the perfect honeycomb billiard, for small and large disorder. 
The symbols are results from simulations, the dashed line is 
a fit, and the full line is based on an analytic result for the 
diffusion constant [see Eq. 1241 1. The channel width is d — 0.1. 
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FIG. 5: Time exponent 7(g) for the gth moment of the par- 
ticle distribution function plotted versus q. The symbols are 
results from simulations for the exact honeycomb billiard with 
channel width d — 1. The fitted dashed lines have slopes 
1.68/2 and 1.72/2. 
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FIG. 4: Particle distribution for t = 10 s for (a) the hon- 
eycomb billiard and (b) the random channel billiard with 
Sr = 0.3. The channel width is d = 0.1. 



in section Hvl For small disorder (Sr = 1CP 6 , star sym- 
bols), the mean-square displacement exhibits a transient 
superdiffusive regime for small times with the same ex- 
ponent v = 1.72 as in the regular case and then enters 
conventional diffusion (y = 1) for large times. Finally, 
for large disorder (Sr = 0.3, circle symbols), the motion 
is purely diffusive. 

Figure 0Jl) shows clearly that the superdiffusive mo- 
tion in the honeycomb billiard is associated with a non- 
isotropic probability distribution P(r, t) of the particles. 
The spikes in P(r,t), plotted for t — 10 5 , suggest that 
the long effective steps, responsible for superdiffusion, 
occur along the six equivalent directions of the channels. 
Within the theory of continuous-time random walks, we 
can show that such a spiky shape of the distribution has 
to appear. However, in the regime of conventional dif- 
fusion (Sr — 0.3), the distribution P(r,t) assumes the 
expected isotropic shape of the Gaussian distribution, as 
illustrated in Fig. 0Jd). 

We investigated if the probability distribution P(r, t) 



P(r,t) = —p(*l 



(1) 



If it is valid, the moments of P(r, t) fulfill 

(|r(t)|«) = (|r(l)H with 7(9) - uq/2 , (2) 

as one can show in a straightforward manner. We deter- 
mined the exponents "f(q) from a double-logarithmic plot 
of (|r(i)| 9 )/(|r(l)| <? ) versus time t. Figure |5] plots -f(q) as 
a function of q for the honeycomb billiard for a chan- 
nel width d — 1. Since the single points follow a nearly 
straight line, the scaling law of Eq. is roughly fulfilled. 
A closer inspection, however, reveals that the regions for 
q < 2 and q > 2 are better fitted by different slopes 
v = 1.68 and v = 1.9, respectively. So the exponent 
v = 1.72 determined from the mean-square displacement 
lies between these two values. Interestingly, such a small 
difference of the slopes was also found by Sanders and 
Larralde in Ref. [l(| for their "parallel zigzag billiard" 
with the kink at q — 3. On the other hand, in the infinite 
horizon billiards with (r 2 ) oc that studied by Armstead et 
al. [23j , the same analysis also reveals a kink at q = 2 but 
with a larger difference of the slopes v — 0.5 for q < 2 
and v = 1 for q > 2, respectively. Finally, we note that 
for random channel billiards 7(g) = q/2 confirming the 
expected Gaussian distribution. 

In Fig.[S]we plot (r 2 ) versus time for different strengths 
of disorder. At small times, there is always superdiffu- 
sional behavior with the same time exponent. Both the 
crossover time to conventional diffusion and the diffusion 
constant D increase with decreasing Sr. Especially, we 
find that D diverges for Sr — > 0. 
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FIG. 6: Mean-square displacement as a function of time t 
for different strengths Sr of disorder illustrating the transient 
regime of superdiffusion. The full and dashed lines indicate 
the limiting cases of conventional diffusion and pure superdif- 
fusion. The symbols are results from simulations, the channel 
width d is 0.1. The inset shows the mean-square displace- 
ment calculated with the step-time distribution of Eq. 11811 
for different tdis- 



IV. SUPERDIFFUSION: AN ANALYTIC 
APPROACH 

A. Exact Honeycomb Billiard 

In Fig.[7|we plot randomly chosen particle trajectories 
from our simulations. They demonstrate that particles 
move, on average, in one direction for a very long time 
before they change their course. Moreover, long steps 
along one of the six possible channel directions are in 
the majority. These long steps occur due to the exis- 
tence of ballistic trajectories or perfect paths, as we call 
them, where particles follow forever a certain direction. 
Examples were already introduced in Fig. [21 The main 
contribution to superdiffusion comes from paths that are 
close to perfect, i e., paths whose initial conditions differ 
slightly from a perfect one. Almost perfect and perfect 
trajectories then take the same channels for a long time 
until the difference between them is so large that they en- 
ter different channels. Therefore, an almost perfect path 
follows the direction of a ballistic trajectory for a long 
time before changing its direction. In the following, we 
consider the straight sections of an almost perfect path 
as an effective step in a Levy walk and use the formalism 
of a continous-time random walk to determine the mean- 
square displacement and the particle distribution func- 
tion on the basis of a distribution for these effective steps 
[l9|. Since in our case, the long steps are restricted to 
special directions, the exponent in the mean-square dis- 
placement is equivalent to the one found in Ref. |l7j for 
pure one-dimensional rather than two-dimensional sys- 
tems, as we will demonstrate explicitely below. 

In a periodic Lorentz gas with infinite horizons the situ- 
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FIG. 7: Simulated paths of particles started with random 
conditions in a perfect honeycomb lattice. Lengths are given 
in units of the edge length io of a hexagon. It is obvious that 
some particles follow straight paths along the six lattice axis 
for very long times before they change direction. 

ation is similar. The main contribution to superdiffusion 
is due to paths close to the straight ballistic trajectories 
depicted in Fig. QJi) . An effective step in the Lorentz gas 
ends when a particle whose initial conditions differ from 
the ones of a ballistic path hits a scatterer. For both the 
Lorentz gas and the hexagonal billard, let us denote by 
e the difference in the starting angle of a perfect and an 
almost perfect path. Particles starting at the same point 
will accumulate a difference Sx in position that in leading 
order grows as te, where t is the travel time of the parti- 
cles. An effective step ends if Sx exceeds some treshold 
Sx m ax proportional to the distance of scatterers in the 
Lorentz gas or to the channel width d in the honeycomb 
billard. In both systems, the duration t of an effective 
step is therefore proportional to 1/e for sufficiently small 
e. For the Lorentz gas this relation follows from the work 
in Ref. 0; for the honeycomb billiard we explicitely jus- 
tify it for a special class of perfect paths in Appendix 

El 

Let p (e) de be the probability of finding a particle close 
to a perfect path with a starting angle in the intervall 
[a + e, a + e + de], where a is the starting angle of the 
perfect path. With t oc 1/e, one finds the distribution 
function A(i) of the effective step times: 

X(t)dt = p (e) de oc ^^-dt. (3) 

In a Lorentz gasp (e) oc sin e, i. e., the velocity component 
cos e perpendicular to the direction of the infinite hori- 
zon is equally distributed Therefore, X(t) oc l/t 3 in 
the long-time limit that leads to a mean-square displace- 
ment (r 2 ) proportional to tint H, 0| characterizing a 
marginally anomalous diffusion. 

In our simulations of the particle transport in the hon- 
eycomb billard, we choose all starting angles with the 
same probability. A constant angular distribution p{e) 
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leads to X(t) oc l/t 2 and therefore (r 2 ) oc t 2 /\nt [l^ . 
However, in the simulations we find ^r 2 ) oc t 172 . We 
suspect that during the evolution of the photon cloud 
the angles close to perfect paths are no longer equally 
distributed, as also observed in the periodic Lorentz gas. 
In the following, we therefore consider an angular distri- 
bution p(e) oc e' 3 with an exponent (3 between and 1 in 
Eq. |PJ . To obtain a normalizable step-time distribution 
X(t), we use a small cutoff at a time r and arrive at 



X(t) oc 



t 2+f3 



'(t-r) 



(4) 



where the step function 6 (t — r) is one for t > r and zero 
for t < t. In the formalism of continuous-time random 
walks, the important quantity is the probability ip(r,t) 
that the particle performs a step along r in time t. It is 
written as ^(r, t) = p(r|£)A(i), where p{v\t) denotes the 
conditional probability to move along r in time t. As sug- 
gested by Fig. [7] we assume that the particles move with 
an effective velocity v along the six possible channel direc- 
tions given by the unit vectors e.,- = (cos7rj/3, sin7rj'/3). 
So the conditional probability becomes 



1 6 



(5) 



We now calculate the particle distribution function 
P(r, t) and the mean-square displacement using the "ve- 
locity model" for Levy- Walks introduced in [ljj . Within 
the formalism of continuous-time random walks, the par- 
ticle distribution function, i. e., the probability to find a 
particle at location r and time t, is given by the integral 
equation: 



P(r,t) = d 2 r' dt'P(r',t')^(r-r',t-t') 



+R(r,t), 



(6) 



where i?(r, t) is the probability to reach or pass the point 
r at time t within one step: 



R(r,t)=p(r\t) / dt'X(t'). 



(7) 



Performing a Fourier transformation in space and a 
Laplace transformation in time, Eq. © can be solved: 



P(k,u) 



R (k, u) 
l-?(k,u)' 



(8) 



Here the bar indicates a function in Fourier-Laplace 
space. Using ip(r,t) = p(r\t)X(t) together with the re- 
spective definitions J5J and ijjj forp(r|i) and R(r,t), we 
find: 



_ X)j=i [l — A (m — ivk ■ ejj\ /[u- ivk ■ e^] 

I { K!, XL i — t. ~ zz " . 

Ej=i[l-A(t»-<«k- ei )] 



(9) 



The Laplace transform of the step-time distribution Q 
can be expanded for small u; 



X(u) 



bu, 



(10) 



where a = V (—1 — /3) and b = (1+(3)t/(3 are positiv con- 
stants and r(x) denotes the Gamma function. Note that 
the normalization X(t)dt — 1 is fulfilled by A (0) = 1. 
Then, for small u and k, Eq. JUJl becomes 



P(k,u) 



T 6 


—a (it — zi>k • ej)' 3 + 6 






—a (it — ii>k • e^) 1 " 1 ^ + &u 



(11) 



In the denominator, the linear term in k has vanished 
because of Y^=o k ' e i = 0- 

The Laplace transform of the mean-square displace- 
ment can be calculated with the help of P (k, u) : 



-V 2 P(k, M )| k=0 . 



(12) 



We are therefore taking a closer look at Eq. J3J in 
the limit k/u — > 0. We expand (it — iuk • ej)' 3 and 
(it — it>k ■ ej) 1+ ^ in terms of k/u and finally obtain 



P(k,„)« 1-1^^/5-3 

it 2 o 



(13) 



Using this approximate form in Eq. i|12|) and calculating 
the inverse Laplace transform, we find the mean-square 
displacement in the long-time limit 



(r 2 (t)) oct 2 -' 3 . 



(14) 



In the region < j3 < 1, to which our calculations ap- 
ply, the time exponent v — 2 — (3 varies between one 
and two so that we are in the superdiffusive but subbal- 
listic regime. As already mentioned above, the relation 
between the exponents in the step-time distribution X{t) 
and the mean-square displacement is the same as the one 
found in Refs. 0, 0, for one-dimensional systems. 
The exponent v = 1.72, which we observe in our simula- 
tions, is achieved with an exponent (3 — 0.28 in the an- 
gular distribution p(e) oc e' 3 . Our numerical results can 
therefore be explained with the assumption that close 
to perfect paths the density of possible particle paths 
sharply drops to zero. 

Fig.|S]shows P (k, u) as given in Eq. 10 for u = 10~ 3 , 
v = 1 and (3 = 0.28. One clearly sees a six-fold symmetry 
that through the inverse Fourier-Laplace transformation 
is also visible in P(r,t). To analyze this star-like dis- 
tribution pattern further, we investigate P (k, u) in the 
limit u/k — > 0. Expanding Eq. Q in u/k and employing 
polar coordinates, k = (k cos tpk, ksmtpk), we obtain 



P(k,u) 



U u 



(15) 



where $ (ipk) is a function that only depends on the an- 
gular variable: 



*(Vfc) 



1 6 7T 

6^[ iCOS (^~f J ) 



1+0 



(16) 
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FIG. 9: Angular part $ (<pk) of the particle distribution func- 
tion in Fourier-Laplace space in the limit u/k — > for /3 = 0.28 
[see Eq. fTM 



FIG. 8: Fourier-Laplace transform of the particle distribution 
function, P(k, u), as given by Eq. for u — 10~ 3 , M = 1 
and P = 0.28. 



Due to the six- fold symmery, $ (^fe) is a real function 
and can also be written as 



1 



.5 



1=1 



1+/3 



(17) 



We plot in Fig. El for f} = 0.28. It has minima 

at cj)k = J7t/3 (j = 1...6), indicating that along these 
directions the width of P (k, it) is smaller compared to 
other directions, as is also visible in Fig. El Therefore, in 
real space the width of the particle distribution will be 
largest along the corresponding six channel directions. 
This explains the non-isotropic distribution pattern we 
find in our simulations [see Fig. 0^.)]. So the star-like 
structure is indeed a result of the fact that effective long 
steps only occur along special directions. 



B. Random Channel Billiard 

In a disordered channel system, perfect paths do not 
exist. However, for small disorder, there are particle tra- 
jectories that pass the same channels as a perfect path 
until they ultimately take another route. This behavior 
is similar to the almost perfect paths in the honeycomb 
billiard. Small disorder therefore cuts off very long steps. 
To take this effect into account, we introduce an expo- 
nential factor in the step-time distribution: 



considerably reduced in number. To arrive at the correct 
limit of the ordered honeycomb billiard, tdis — * oo for 
5r — ► 0. On the other hand, our description based on 
the step-time distribution (|18|) breaks down, when tdis 
approaches one, i.e, the time a particle needs to pass one 
channel. Diffusion in such strongly disordered channel 
billiards can be modeled by assuming that, on average, 
the particle's paths in different channels are uncorrelated, 
as we will explain below in detail. 

The new step-time distribution (|18|l has a finite first 
and second moment. Therfore, the mean-square displace- 
ment in a disordered channel system is linear in time in 
the long-time limit, i.e., the particles behave diffusively. 
For smaller times, superdiffusion occurs depending on the 
value of tdis, as illustrated by the inset of Fig. At a 
time of the order of tdis, a transition to the diffusive 
regime occurs. This is in aggreement with the results 
from our simulations. For decreasing disorder, the cross- 
over time tdis becomes larger and so does the time range 
where superdiffusion is found. Ultimately, in the limit 
Sr — > the time tdis diverges and superdiffusion persists 
in the long-time limit. 

At sufficiently strong disorder in the random channel 
billiards, the particle's paths in different channels are not 
correlated on average. Then we can consider a random 
walker that performs steps along the edges of the Voronoi 
tesselation. On average, the step length is equal to the 
edge length of the unperturbed honeycomb lattice, i.e., 
one in our scaling [27^ . Furthermore, after each step the 
random walker changes its direction with the same prob- 
ability 1/2 either to the right or to the left by an angle 
that, on average, assumes a value of 7t/3. This enables 
us to calculate the mean-square displacement: 



e -t/t dis 

\ dis (t) cx e-'/* d "A(i) cx t2+fj {t-r), (18) 

where tdi S is a time that decreases with increasing disor- 
der dr. Paths with effective step times beyond tdis are 



( r2 ) = (E r ^^) = ?i + 2 EE^-^ 



= 1 3 = 1 



(19) 



where the vector characterizes a single step with |r$ 
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1. It is straightforward to show that (r^) = Yj/2 l ~i for 
i > j and we can set (r^ • r 3 -) = r 3 • r 3 /2 1_J = l/2 i_J . 
The mean-square displacement is then calculated using 
the formula for geometric sums: 



\ / / j / i Oi—j On 



In the limit of : 



t=i 3=1 

■+ oo this becomes 
(r 2 ) = 3n. 



(20) 



(21) 



To arrive at the mean-square displacement in terms of 
time t, we have to relate the number n of steps to t. We 
first consider a particle path with a fixed angle a relative 
to the channel direction. For small channel widths d so 
that details at the channel junctions can be neglected, 
the time T(a) to pass the channel is 



T(a) 



1 



cos(a) 



(22) 



The step number n is then given by the average over all 
possible angles a: 



2 r /2 t , 2 

71 = ~ / ^FT^ da = ~t ■ 



(23) 



With Eq. (|2~TJ) and (r 2 ) = 4Dt, we finally arrive at the 
diffusion constant 



D = 



2tt' 



(24) 



which is an excellent estimate for random channel bil- 
liards as illustrated by the perfect fit of the full line in 
Fig-Elto the simulation results for 6r — 0.3 and d = 0.1. 
Note that Eq. (|24|) does not depend on disorder. As soon 
as all long-time correlations are destroyed, any further in- 
crease of disorder does not affect the diffusion constant. 



V. CONCLUSIONS 

With the honeycomb billiard that can also be viewed as 
a system where hexagonal scatterers are placed on a tri- 
angular lattice, we have investigated a periodic extended 
billiard in detail. Though particles moving in the hon- 
eycomb billiard always have a finite horizon, there exist 
perfect paths where they move ballistically in one direc- 
tion. We have clarified that almost perfect paths give rise 
to an overall superdiffusive behavior. Our simulations re- 
veal a mean square displacement (r 2 ^ cx t v with a time 
exponent v = 1.72. On the other hand, in our analytical 
treatment we have applied the Levy walk model based on 
the formalism of continuous-time random walks by con- 
sidering the long straight parts of almost perfect paths 
as effective steps. Assuming for their occurence a gen- 
eral distribution of the form p(e) cx e' 3 , we can show that 



the mean square displacement possesses a time exponent 
v = 2 — /3. A comparison with the simulation results then 
gives (3 — 0.28 which means that almost perfect paths are 
less probable than other trajectories. 

In contrast to previous treatments of Levy walks, the 
directional distribution of steps in our model is not 
isotropic. Instead, steps are limited to the six possible 
channel directions, which are used by most of the perfect 
paths as suggested by our simulations. We therefore find 
that the time exponent of the mean square displacement 
corresponds to the one determined for one-dimensional 
systems ^3,0>E|- Furthermore, our analysis reveals a 
starlike distribution of the particles' positions with a six- 
fold symmetry in accordance with our simulation results. 
So a limitation of the allowed step directions together 
with a step-time distribution that causes superdiffusion 
ultimately gives rise to an anisotropic particle distribu- 
tion. 

We have also introduced disorder into the honeycomb 
billiard so that the directions of the channels are ran- 
domized. In the limit of long times, these random chan- 
nel billiards always display diffusive behavior. Transient 
superdiffusion is, however, visible in systems with weak 
disorder for small times. We explain it with the help 
of an exponential cut-off in the step-time distribution. 
For large disorder, correlations along the particle path 
between successive channels are lost on average. With 
the help of an elementary random walk model on a hon- 
eycomb lattice, we can estimate the diffusion constant 
which fits our simulation results very well. 

Our investigation shows that the transport of particles 
in two-dimensional channel systems governed by the rule 
of specular reflection adds further insight to the current 
knowledge of extended billiards. Since we were moti- 
vated to the present study by light transport in foams, 
as mentioned in the introduction, the interesting ques- 
tion arises how the superdiffusion of classical particles in 
the honeycomb billiard will affect the analogous problem 
of interfering waves travelling along the channels. This 
question is also crucial for the relation between classical 
and quantum mechanics psj . 
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a 



i+i 



case b 



0<cc<7t/6 




=71/3-01; 
ji/6<a j + [<Ji/3 



a i+2 =|it/3-a i+1 |=a ; 
0<a ; , 9 <7c/6 




FIG. 11: All possible angles for a particle with a starting 
angle on < 7r/6. In every second channel the starting angle is 
repeated: a i+ 2 n = a,. 



FIG. 10: The position of the particle within the ith channel 
is given by the distance \xij\ or \xi.i\ between vertices of 
the underlying honeycomb structure (dotted lines) and the 
first or last intersection of the particle's path with the center 
line of the channel. The sign of Xi t i and Xij is negativ if the 
intersection is to the left of the vertices otherwise it is positive. 
The angle of the particle's path with the center line is denoted 
Qj. It is always taken positive and does never change within a 
channel. For a,: < 7r/6, the particle proceeds into the channel 
i + 1 situated either opposite (case a) or next (case b) to the 
last reflection in channel i. 



APPENDIX A: DISTRIBUTION OF EFFECTIVE 
STEP TIMES 

In this appendix we calculate the step-time distribu- 
tion for our Levy-walk model in Sec. IIV Al In concrete, 
we will consider paths with effective steps along the main 
lattice directions. 

For special starting angles a%, so-called perfect paths 
exist where the particles run forever along one of the 
lattice illustrated in Fig. [21 A particle having 

started at an angle ai + e will leave such a perfect path 
after an effective step time t. Analyzing special perfect 
paths to be defined below will allow us to calculate the 
step-time distribution X(t). 

We first take a closer look at a general path of a particle 
crossing the zth channel. The position and direction of 
the particle is characterized by Xij (or x^i) and ai, as 
defined in Fig. EH Depending on the position, a path 
with an angle a, < 7r/6 proceeds into the new channel 
i + 1 situated either opposite (case a in Fig. ITU|) or next 
(case b) to the last reflection in channel i. The new, 
respective angles therefore are a^+i = 7r/3 — a, or a,+i = 
7r/3 + ai (see Fig. Illf) . However, in channel i + 2 one 
always finds ai+2 — ca since for aj+i > tt/6 only case a 
applies. 

For reasons to become clear below, we now calculate 
the position 2Ji+2,/ in channel i + 2 as a function of the 
parameters a, and x% f in the ith channel assuming a^ < 
7t/3. We know already that this is fulfilled in every second 
channel. First, we calculate the position x^i from the 
position Xij: 



where the integer mi is the number of reflections in chan- 
nel i. Secondly, with the law of sines we determine x i+ ij 
(see Fig. EH: 



sin ai 

x i-\-l,f — ~F — : -EiAast- 

sinQi + i 



(A2) 



The upper and lower sign belong, respectively, to case a 
or b. The relations for channel i + 1 equivalent to Eqs. 
(|AlT) and HA2j) are 



Xi+i,i = Xi+ij + m i+ idcot Oj+i - 1, (A3) 

(A4) 



sin a l+1 

Xi+2J = T ZIZ I x i+l,l> 



sin a.i 



where in Eq. (| A4|> a;+2 = ai was used. Finally, combin- 
ing all equations IjAlfl to i|A4(l . we obtain: 



%i+2,f = Xij — 1 + Widcot ai 

sin (? qp a i) cos(^=Fai 



±- 



sm ai 



T 



sin at 



■m l+ id.(A5) 



From the multitude of possible perfect paths, we cal- 
culate the step-time distribution X(t) for a special class 
of perfect paths characterized by the requirement that 
the position Xij in the channel is periodically repeated 
in every second channel, i. e., Xi+2kj = %ij for all k. 
From Eq. I|A5(1 . this is the case for 

tana, tana 4+ i 

mi = ; — and m i+1 = , (A6) 

a a 

where now belongs to this perfect path. Examples are 
illustrated in Fig. E with paths (1) and (3). 

For a particle with starting angle a^ + e, not running on 
a perfect path, Eq. I|A5J| applies as well with ai replaced 
by ai + e. We can therefore calculate the difference in 
positions 8xi + 2 — ^+2 / — x i+2,f m channel i + 2 when 
the difference in channel i is Sxi = x^ * — x^j: 



8x i+2 = Sxi + K ai 



sin (ai + e) ' 



where 



= Xij + rrii d cot ai 



1 



(Al) 



K„. = 



=F3 — v^tanai 
± cos ai + v3 sin ( 



(A7) 



(A8) 
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If Sx n exceeds a threshold proportional to Sx max = 
dj tanai, particles travelling on the perfect path and its 
neighbor proceed into different channels (see Fig. If 2JI . 
Therefore the maximum number of channels n max where 
they travel through the same channels is 



sin («! + e) / tanai sinai ' 

oc ox max : = d smai H 

sm e \ tan e 

(A10) 

Since the particles move ballistically and for small e, the 
duration of an effective step therefore is 



'max 



FIG. 12: Two particles can only enter the same channel if the 
distance between the last intersections of their paths with the 
channel center line is smaller than 5x ma x = d/ tana n . 

and Eqs. ljA6() were used. We now assume that particles 
on the perfect path and its neighboring path start at the 
same position, i. e., 8x\ = 0. In the long-time limit the 
different behavior in channel i and i + 1 is irrelevant so 
that in the nth channel the difference in position Sx n 
becomes 

ox n oc n ; -. (A9) 

sm (ai + e) 



tocn max (x-. (All) 
e 

The number of particles with starting angles in the in- 
tcrvall [ai + e,ai + e + de] is p(e)de, where p{e) is the 
distribution of starting angles close to a perfect path. 
Together with Eq. (| Al 1|) . this determines the step-time 
distribution for t — > oo: 



X(t)dt = p (e) de oc , y dt. (A12) 
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